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ABSTRACT 



Cassini radio science experiments have provided multiple occultation optical depth profiles of 
Saturn's rings that can be used in combination to analyze density waves. This paper establishes an 
accurate procedure of inversion of the wave profiles to reconstruct the wave kinematic parameters 
as a function of semi-major axis, in the nonlinear regime. This procedure is established using 
simulated data in the presence of realistic noise perturbations, to control the reconstruction error. 
It is then applied to the Mimas 5:3 density wave. 

There are two important concepts at the basis of this procedure. The first one is that it uses the 
nonlinear representation of density waves, and the second one is that it relies on a combination of 
optical dept h profiles instead of just one p rofile. A related method to analyze density waves was 
devised by lLongaretti and Borderies (1986) to study the nonlinear density wave associated with the 
Mimas 5:3 resonance, but the single photopolari metric profile provided limited constraints. Other 
studies of density waves analyzing Cassini data (|Colwell and Espositoll2007 : Tiscareno et al. I l2007h 
are based on the linear theory and find inconsistent results from profile to profile. Multiple cuts 
of the rings are helpful in a fundamental way to ensure the accuracy of the procedure by forcing 
consistency among the various optical depth profiles. 

By way of illustration we have applied our procedure to the Mimas 5:3 density wave. We 
were able to recover precisely the kinematic parameters from the radio experiment occultation 
data in most of the propagation region; a preliminary analysis of the pressure-corrected dispersion 
allowed us to determine new but still uncertain values for the opacity (K ~ 0.02 cm 2 /g) and 
velocity dispersion of (c ~ 0.6 cm/s) in the wave region. 

Our procedure constitutes the first step in our planned analysis of the density waves of Saturn's 
rings. It is very accurate and efficient in the far-wave region. However, improvements are required 
within the first wavelength. The ways in which this method can be used to establish diagnostics of 
ring physics are outlined. 



Key Words: Saturn's rings, nonlinear density waves 
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1 Introduction 



The wealth of new data provided by the Cassini mission (and in particular, the collection of in- 
dependent radial wave optical depth profiles of Saturn's rings) opens the possibility to use wave 
dynamics to establish physical diagnostics of ring physics to a level of accuracy substantially higher 
than what was possible with Voyager. However, such a research program requires an accurate de- 
termination of the wave kinematic parameters throughout the wave propagation region. This has 
prompted us to reinvestigate this last problem. 

Many density waves in Saturn's rings are strong and nonlinear - tha t is the y cannot be mod- 
eled by the linear theory of density waves fully described e.g. by IShul (|1984|) - but this linear 
theory is nevertheless used due t o the absence of a w ell-established non l inear inversion proce- 
dure (e.g., iNicholson etaDll990l : iRosen et all 1 1991 al ibi: ISpilker et al]l2004l : iTiscareno et al.ll2006l 
20071) . This paper presents a new approach to analyze linear and nonlinear density waves. The 



procedure is established through the use of a number of simulated optical depth profiles, as only 
simulated data allow us to compare the reconstructed kinematic parameters with the original ones. 
It is applied next to the Mimas 5:3 density wave as observed by the Cassini Radio Science oc- 
cultation experiment, as an example. This procedure will serve as the basis of several subsequent 
studies of real data. The approa ch adopted here considerably sharpens the method described by 
Longaretti and Borderiesl (|l986l) . which extracted the maximum information from one Voyager 



photopolarimeter stellar occultation profile. However, one profile is insufficient to accurately con- 
strain all the parameters of the model to the level of accuracy required to produce new and detailed 
diagnostics of the ring physics. 

There are several motivations for analyzing density waves, and especially nonlinear density 
waves. The first objective of determining the kinematic behavior of ring density waves is to test 
the background kinematic model described in section[2l as possible systematic deviations from this 
model might provide information on physics that are yet unknown or not well constrained. Cassini 
provides many different radio occultation profiles of the rings, which can be used to verify that 
a single set of wave parameters can be determined for various observations of a wave at various 
longitudes with respect to that of the associated satellite (at least when modulations due to satellite 
orbital variations or other physical effects are weak enough). 

For each wave, the set of kinematic parameters determined by the method described in this 
paper, coupled to the analysis of the evanescent part of the wave (which may require extrane- 
ous dynamical constraints to be accomplished), form the basis of the next objective, which is to 
measure as precisely as possible the torqu es exerted by the satellites on th e rings. Such a direct 
measurement was previously attempted in lLongaretti and Borderiesl (| 1986Q . but with limited suc- 
cess due to insufficient constraints in the first wavelength of the studied wave, and to failure of 
the WKBJ approximation in the vicinity of the resonance. This type of torque measurement based 
on a determination of the wave kinematics can only be performed in Saturn's rings, and therefore 
constitutes a unique way to verify the dynamical theoretical predictions for such torques. The 
measurement and verification of the torque are directly relevant to the dynamics of disk-satellite 
interactions in general. The extent to which dynamical constraints are required on top of a purely 
kinematic description of the wave remains to be seen; this bears on the model dependence of such 
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a torque measurement. 

Another important study that will follow from this paper is the detailed analysis of the ring 
stress tensor, with the hope that this may in turn provide useful constraints on the ring particles 
collisional properties. We wish to identify statistically the various stress behaviors that might 
be found, depending, for instance on th e ring region or ring b ackgr ound optical depth . Several 
models of damping exist. For instance, iBorderies et al. I dl983b and IShu et all dl985ah predict a 



bimodal behavior as a function of the ring mean optical depth in dilute rings. IBorderies et al. 



dl985h generated models for dense rings. This research program requires combining the kinematic 
reconstruction method developed here with the nonlinear (pressure-corrected) dispersion relation 
and at least one ge neric dynamical e quation describing wav e propagation and damping (such as 
the ones derived in lshu et al]|l985al or IBorderies et al.lll986|) . 

In the bulk of this paper, we simulate a density wave that has some similarities with the Mimas 
5:3 density wave, as observed in eight diffraction-corrected Cassini radio occultation profiles. We 
did not try to simulate this wave precisely, because the purpose of the simulation was to establish 
the procedure. A preliminary and simplified analysis of real data pertaining to this wave is pre- 
sented in section \5\ The reader might wonder why we started with simulated data. The reason is 
that such an analysis is in itself already very complex and it is important to separate this complexity 
from that inherent to the analysis of real data. There are several reasons for which the analysis of 
simulated data is difficult. 

The first reason is the nonlinearity in itself. The equation describing the shape of the wave, 
with its large troughs and narrow peaks, does not provide a model that one can readily fit to the 
data. The fit turns out to be extremely sensitive to the dependence of the wave phase on semi-major 
axis, and so the procedure must allow for a very precise determination of this quantity. By using 
simulated data, we can compare how the reconstructed wave kinematic parameters such as the 
wave phase compare to the ones used in constructing the data; we can therefore quantify the error 
in the reconstruction method developed here. 

The next difficulty is associat ed with the fact that the absolute radial scale of Saturn's rings i s 
presently known to only ±2 km dNicholson et al.lll990l : French et al.lll993c IJacobson et al.l l2006). 
Long period variations in the satellite's orbit may produce similar radial shifts in the resonance lo- 
cation, and are expected to produce variations in the wave propagation which are not yet modeled. 
The uncertainty in the radial scale, at least, will be reduced by future kinematic models to the level 
of about 100 m. 

Another major problem comes from the fact that one must distinguish between radius r and 
semi-major axis a; the relation r = a [1 — ecos(m0 + mA)] (see section [2] for definitions of the 
variables and discussion of this relation) is simple only superficially. 

Other effects add their own complexity to the data. Interestingly, the noise present in real data 
(which we took into account in our simulations) is not a real problem, although the need for several 
profiles to constrain a unique self-consis tent solution result s mostly from the presence of noise. 
However, and as indicated by the work of iLewis and Stewart! (|2000, 2005), gravitational clustering 
can cause significant effects on the structu re of local density maxima. Also, the gravita tional 
wakes observed in the rings (see for instance Hedman et al.ll2007l : IColwell et al.ll2006Ll2007|) make 
it necessary to normalize the background optical depth among the various profiles, and may affect 
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the optical depth profile of the wave in an as yet uncontrolled way. Focusing on simulated data 
allows us to ignore this added complexity in the first step . 

Before describing our procedure in section H we begin in section [2] by gathering the basic 
equations of the theoretical streamline kinematic model. Section [3] describes and illustrates the 
method to obtain the simulated data. Section [5] applies the procedure to the Mimas 5:3 density 
wave. Section [6] summarizes and concludes the paper with a discussion of future work. 



2 Theoretical Model 

The most striking observed characteristic of density waves - the shape of the peaks and troughs - 
is kinematic in nature. Dynamics enters mostly through the nonlinear dispersion relation, which 
specifies the change of the wavelength with radius as the wave propagates, and through the control 
of the wave amplitude; another dynamical feature is the existence of an evanescent zone. More 
generally, dynamics controls the way the observed kinematic features are spatially modulated. 
It is therefore legitimate first to focus on the information that can be recovered from kinematic 
constraints, relying as little as possible on dynamical ones. 

The major difficulty of this program lies in the fact that the observed wave profiles depend in a 
complex, mostly implicit way on the underlying kinematics, making the elaboration of an inversion 
procedure both a necessary and complex task. In this section, we gather the required kinematics 
(and the minimal dynamics) needed to fulfill this objective. 

This section addresses test particle kinematics (section [27Tb . streamlines (section |272~1) . surface 
density and optical depth (section 12.31) . and finally discusses the link between radius and semi- 
major axis ( section |2~4l ), which plays an important role in the inversion procedure. Particular atten- 
tion is paid to the discussion of the generality of the kinematic model and dynamical constraints 
used in our inversion procedure. 



2.1 Test Particle Kinematics 



We consider a satellite orbiting in the equatorial plane of the planet. The cylindrical coordinates are 
denoted by (r, 6), where r is the radius and 9 is the longitude. We use the subscript s to characterize 

the satellite; 

As in lGoldreich and Tremaind (|1982i ). one can linearize the equations of motion of a test parti- 
cle orbiting the planet to obtain 



d 2 r fd9\ 2 <9$ 

dt V dt 0r' 1 J 



where 
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$ = % + $ s (3) 

is the sum of the planet potential and the satellite potential, to obtain the first order solution for r 
andtf. 

At any given time t, the satellite potential exerted on a particle at (r, 6) is periodic in (9 — 6 S ). 
Furthermore, the satellite makes an epicyclic oscillation at angular frequency k s , the epicyclic 
frequency, about a guiding center which revolves at the rate Q s , the mean motion. Thus the satellite 
potential can be expanded in a double Fourier series, one in longitude and one in time: 



-co + OO 



$ s (r, 9,t) = J2 Yl cos H# - Os) + kK a (t - t )] , (4) 

m=0 fc=— oo 

where the integers m and k characterize each term of the series and to is the time of satellite 
periapse. The pattern speed 

k 

f2 p = Vl s -\ k s , (5) 
m 

is the angular speed of the rotating frame in which the (m, k) component of the satellite potential 
is stationary. 

Denoting by Vto the angular rate of the ring particle in the zero order solution of Eqs. (OQ) 
and © for r and 6, the first order solution is singular if Q p = tt (corotation resonance) or if 
m(fi — tip) = ±k (Lindblad resonance). Inner Lindblad resonances correspond to the + sign and 
occur in a ring inside the satellite, and outer Lindblad resonances correspond to the — sign and 
occur in a ring outside the satellite . We will restrict ourselves here to inner Lindblad resonances. 
The resonance is labeled (m + k) : (m — 1) because we have 

Q m + k 

(6) 



Q q m — 1 



an approximation that we note here but will not use in the remainder of this paper. 
Near a Lindblad resonance, the linear response of a test particle is given by: 



r = r + -7— cos mcj), (7) 
Ar 

9 = 9 + tt t + -^sinm0, (8) 
Ar 

where Ar is the distance from the resonance and 

ni^ = m [6 - 6 s0 - Q p (t - t )} , (9) 

with 

0so = Os(t o ). (10) 
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The coefficients A r and A e are given explicitly in lGoldreich and Tremaind (I1982f) . 

As noted earlier, the amplitude of the first order solution for r and 9 becomes infinite at the 
resonance location. The linearization assumption that the first order solution is much smaller than 
the zeroth order solution breaks down too close to the resonance. 



2.2 Streamlines 

When considering fluid particle motions instead of test particle ones, it turns out that collective 
interactions (collisions and self-gravity) prevent the divergence just mentioned from occurring; 
they also make the description in terms of first order deviations from circularity valid to a high 
level of precision (relative deviations are typically of the order of 1CT 5 or 1CT 4 in observed density 
waves, an order of magnitude that can be deduced either from linear theory or from order of 
magnitude analysis of the data; see lLongaretti and Borderies 19861 for details), while changing the 



amplitude and phase of this first order description. Because waves are forced by external satellites 
and because self-gravity acts as a cohesive force on wave motions, fluid particles sharing the same 
semi-major axis are expected to follow the same m-lobe orbit, or streamline, in the frame rotating 
with the pattern speed associated with the considered satellite resonance [see Eq. ©]. The shape 
of ring streamlines is given by: 

r = a [1 — e(a) cos (m</> + mA(a))] , (11) 

where a is a semi-major axis, e(a) « 1 is the eccentricity, and A(a) is a lag angle; as in the 
previous section, (f> is the azimuth in the rotating frame. The lag angle has a direct geometric 
interpretation: it is the angle by which the m-lobe shape of the wave rotates when moving from 
one streamline to the next [see Fig. \T\ where ip is defined by Eq. ([22])] • 

This description implicitly assumes a Lagrangian approach to fluid motions. An unperturbed 
fluid particle follows a circular orbit, and has coordinates a, <p. Once perturbed, this same particle 
follows an m-lobe orbit and has coordinates r(a, 0), 9{a, 4>), where r(a, <fi) is the relation above 
and 

= (f) + 2— esin(m0 + mA). (12) 

K 

In this description, a, are used as Lagrangian labels of individual fluid particles, instead of 
the more customary initial position r , 9 . It is important to note that a, <p are both the unperturbed 
position of the fluid particle, and its Lagrangian labels. Functions of space, such as the ring surface 
density, are either specified in terms of r, 9 or, more often, as a function of a, <p through the change 
of variables defined by the two previous relations. 

Eqs. (fTTI) and (fT2l) can either be viewed as an a priori assumption concerning fluid particles 
kinematics, to be validated t hrough the derivation o f appropriate dynamical solutions of the equa- 



tions, an approach adopted in Borderies et al. (|l98q) : or it can be derived from ab initio analysis of 



these same equations, as in lShu et al.l (| 1985a ). Note however that Eq. (fTTj) . the most important of 
the two, constitutes in practice the most general such relation one can assume: the fluid particle re- 
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sponse must be periodic in mcp because the satellite forcing i^]; and the magnitude of the observed 
deviations from circularity ensures that a sinusoidal periodic response (the first term expansion of 
any periodic response) will accurately represent the data to their current level of precision. 

For our purposes, Equation (TTTI) involves several implicit features. The first one is that the 
streamlines do not cross. Non-interacting streamlines would cross at locations where dr/da = 0; 
indeed, wherever a difference in semi-major axis 5a would not produce a difference in radius 
5r = dr/da5a at a given azimuth, fluid particles with different semi-major axes would occupy the 
same location, leading to streamline crossing. One 

dr 

J— — = 1 — q cos imcj) + mA + 7) , (13) 

da \4> 

where we have neglected the small term e cos(m0 + mA) and with: 

de 

gcos7 = a—, (14) 
da 

dA 

gsin7 = mae——. (15) 

da 

In Eq. (fT3l . the lag angle between streamlines is defined by Eq. (TTTI) and 7 and q are defined by 
Eqs. (fT4l) and ([TBI) . Equation (fT3l) is a mathematical consequence of Eq. (TTTI) and requires no extra 
physical input. Streamline crossing is prevented as long as q < 1. It is known from observations 
that although e is a very small parameter, q is usually a significant fraction of unity. The physical 
meaning of q and 7 will be further illustrated in the section devoted to the surface density, [231 

From a dynamical point of view, as all forces are small compared to the planet attraction, fluid 
particles must follow perturbed epicyclic orbits, very much like test partic les follow perturbed ellip 



tic (or possibly epicyclic) orbits. This flui d particle epicyclic theory (see Borderi es and Longarett 



19871 Longaretti and Borderieslll99l[ and lBorderies-Rappaport and Longarettilll994f) was devised 



to put the streamline formalis m on a more self- consistent footing. For an ab initio construction of 



the streamline formalism, see lLongarettil (| 1992Q . We have: 



r = a [1 — e cos(6> — w)] , (16) 

with: 

w = wo + w(t — to), (17) 

where w = — k is the fluid particle apse precession rate. 

Equations (fTTI) and (fT6l) are compatible if (to lowest order in eccentricity): 

(m — l)#o — tt — m0 s o — mA = 0, (18) 



'Only the lowest terms in eccentricity are important. 

2 The name J stems from the fact that dr/da is the leading order term in eccentricity in the expression of the 



Jacobian of the change of variable from r, 9 to a, 1 
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and: 



m(fi - n p ) = Q - w. (19) 

Equation (fT8l expresses the constraint that all fluid particles with the same semi-major axis 
belong to the same m-lobe streamline. This is a generalization of the condition that an elliptical 
ring is such that all the particles with the same semi-major axis have the same periapsis angle. 

Equation (fT9l) looks like the resonance condition, but it includes all the perturbations exerted on 
the particles (whose dominant contributions is on the precession rate), in line with the general phi 



losophy of osculating elements; as such it applies throughout the wave region (see lBorderies et al. 



19861) . This description is useful as long as all these perturbations (arising from interparticle colli- 



sions and self-gravity, mostly) are much smaller than the effects of the planet, a condition largely 
satisfied in ring systems. As a consequence, the streamlines do not differ much from the orbits of 
test particles viewed in the rotating frame. Note that as long as explicit expressions of the forces 
are not used, the dynamical content of Eq. (TT9~1) is absolutely generic, as the osculating epicyclic 
theory is mathematically equivalent to the generic form of the equations of fluid motion in the limit 
of small deviations from circularity. As such, this condition is simply a kinematic consequence of 
Eq. (PTTI) . Equation (PT21) is valid to the same degree of generality. 

The angle A [defined by Eq. CCD] is introduced because the streamline orientation varies 
smoothly from the inside to the outside of the resonance. Indeed, in Eq. (fT9l) , the secular drift 
of m{p, — Vt p ) — k that would be due to the planet alone can be compe nsated by self-gravi ty (a 



condition that by itself yields back the wave dispersion relation; see, e.g. jBorderies et al.l ll986 and 



Longaret ti 1992). But since self-gravity is a small force, this compensation takes place only if the 
shift of A between streamlines is large enough, i.e. if we are in the condition of the tight- winding 
approximation: 

ma^ » 1, (20) 
da 

This is why density waves are so tightly-wound in rings. In spiral galaxies, the self-gravity 
dominates the force due to the central bulge, radial excursions are much larger (so that a first order 
deviation from circular motions such as the one adopted here is much less appropriate) and the 
spiral arms are much more open. The equatio ns describing linear waves can be recovered the limit 



where q « 1 (ILongaretti and Borderieslll986|) 



Another useful qualitative dynamical feature relates to the WKBJ ordering, which is usually 
assumed in theoretical analyses of density waves. Let us point out the meaning of this approxima- 
tion for our purposes. A wave possesses three radial scales^: the semi-major axis, the wavelength 
I (note that the radial wavenumber k ~ m(dA/da), a near-equality that holds to leading WKBJ 
order), and the scale of variation of the wave amplitude and background^, which we will denote 

In order of magnitude, for a typical wave, a ~ 10 5 km, / ~ 10 km, and £ > 10/ (as waves 



3 A fourth radial scale is obtained from ae. This scale is smaller than any of the other three, and is not required in 
this discussion. 

4 For instance, the lengths over which the amplitude of the wave increases and decreases - the damping scale, or 



the length of the wavetrain. 
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propagate for a few to a few tens of wavelengths). The WKBJ ordering holds when the variation 
of the wave phase function (mA + 7) with radius is much faster than the variation of the wave 
amplitude, i.e., fc£ 1. Eqs. (PT4l) and (fT5l) imply then that gcos7 ~ ae/£ <C gsin7 ~ A;ae and 
that 7 ~ 7r/2 when the WKBJ ordering holds. Also, as all density waves excited at inner Lindblad 
resonances display decreasing wavelength and decreasing amplitude q with increasing radius (at 
least in the far propagation region for the amplitude), and as q ~ kae, the decreasing amplitude is 
related to a decrease of the eccentricity on scale £. 

The WKBJ ordering (k£ 1) is satisfied in Saturn's rings, albeit to a much lesser degree of 
precision than the tight-winding approximation (ka ^> 1). This motivates us to assume that most 
of the radial variation of the wave is due to the dependence of the lag angle A on semi-major axis. 
We nevertheless devise a correction to account for the imperfection of the WKBJ ordering, as a 
precise determination of the phase function i s essential to the robustness of our procedure; 



No te t hat both the linear (e.g ., lShulll984 and references therein) and nonlinear (e.g. JShu et al 
1985a and iBorderies et al.l ll986) density wave dispersion relations imply that k — > at the reso- 
nance radius, so that the WKBJ ordering must fail there and 7 ~ in the evanescent part of the 
wave and close to the resonance. This fact limits our ability to reconstruct the wave kinematics 
inside the first wavelength, a point that will be further discussed later on. 

For future reference, we stress that several different "phases" are used in this paper, and intro- 
duce some terminology to distinguish them. These phases are the lag angle A, the phase function 
defined by: 

/(o)=mA(o)+7(o), (21) 

and the profile phase ip: 

iP = m<p + f. (22) 

The phase function is a characteristic of the wave (i.e. all the profiles have the same phase 
function) and the profile phase characterizes the location of peaks and troughs in 2D as well as in 
each ID wave profile. The profile phase is the one which may be associated to the usual concept 
of wave phase; e.g., the density maxima in the (a, <fi) plane correspond to ip = 2n modulo 2n in the 
WKBJ ordering, while minima correspond to ip = n modulo 2n (see Fig.Q]). The phase function is 
a useful intermediate quantity which gathers the semi-major axis dependence of the profile phase, 
and whose derivative is, by definition, the radial wavenumber. 

Summary: This subsection (as well as the previous one) has introduced and discussed a number 
of elements of wave kinematics and dynamics. In our reconstruction procedure, the only kinematic 
relation introduced here that we use is Eq. (fTTj) and its mathematical consequence Eq. (fT3l) . along 
with Eq. (fT2~l) . We have argued that they are quite general; furthermore, Eq. (fT2l) is used only in a 
limited way, to justify that the difference between 9 and <p can be ignored, a consequence of e <C q. 
The only piece of dynamics used is the WKBJ ordering, i.e., the fact that the radial variation of the 
lag angle mA is much larger than the variation of amplitudes and background e, q, and S (the ring 
unperturbed surface density). As a consequence 7 ~ ir/2 over most of the wave. As discussed, 
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this assumption is expected to fail inside the first wavelength, and our simulated data incorporate 
this feature. 



2.3 Surface Mass Density and Optical Depth 

As the length scales of the observed wave pattern are much larger than the local ring thickness, it 
is legitimate to adopt a two-dimensional description. Conservation of the mass of a ring element in 
the Lagrangian change of variable r, 9 to a, <p implies that dM = £(r, 0)drdO = £n(o)| J\dad(j) = 
T: (a)dad(j) so that the perturbed surface mass density given by: 



where we have directly identified the Jacobian of the change of variable with its dominant term 
in eccentricity given by Eq. (fT3l) , while noting that J > for q < 1; £ is the unperturbed 
surface density, also called the background surface density. This follows because the ring mass 
SM enclosed between these two streamlines is the same between perturbed and unperturbed states 
(SM = T,SrS9 = Y^SaScj)), and because, to leading order in eccentricity, the Jacobian of the change 
of variable from a, <p to r, 9 reduces to J. 

This relation provides us with another meaning for q: at peaks, a = E G / (1 — q) while at troughs, 
a = Yj j (1 + q). Therefore q measures both the density variation due to the wave propagation, and 
its nonlinearity, and 7 measures the relative contribution of the variations in eccentricity and lag 
angle to this density contrast [see Eqs. (fl4l) and (fT5l)l. Only the lag angle contribution to the density 
contrast has been used in Fig. Q] (consistent with the WKBJ ordering); this figure also illustrates 
that the variation of distance of the streamlines is the direct cause of the surface density variations 
of the wave, consistently with Eqs. (|23l and dT3b . 

Application of Eq. (T23l) to the ring normal optical depth r is usually taken for granted, but 
requires some discussion. For simplicity, assume that the ring particle distribution with respect 
to size, spin, shape, etc. can be represented by a discrete index i, and let iVj(r, 9) be the column 
density of particles of type i, rrii their mass, and Si their cross-section. With these definitions, the 
ring surface density E(r, 9) = A^mj while its optical depth r(r, 9) = ^ NiSi. If the distribution 
is constant in time, the formal similarity between these two relations implies that the optical depth 
r will also formally obey Eq. (T23T) . as mass conservation then reduces to number conservation 
and as the wave restoring force (the self-gravity) does not induce particle segregation. In such a 
context, the optical depth is therefore also given by 



The remarks above show also that the same relation will hold to a high degree of precision if 
the size distribution, in particular, is evolving slowly in comparison to the dynamical time-scale 

5 i.e., the optical depth that would be found if the deviations from circular motions were suppressed, but this may 
differ from the optical depth outside the wave propagation region; see the following discussion. 



So 



(23) 



1 — q cos (m0 + mA + 7) ' 



T = 



1 — q cos (m0 + mA + 7) ' 



(24) 
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of the wave. This depends i n turn on the efficiency of gravitational accretion in the ring region. 



Weidenschilling et al.1 41984) put forward and studied the concept of Dynamical Ephemeral Bodies, 



or DEBs. T his concept l eads t o a particle size distribution evolving on the time- scale of a rotation 
period 1 dLongaretti . 1989|) . but requires that the efficiency of gravitational accretion is large 
and only mildly variable throughout a sizeable frac t ion o f the ring to account for the remarkably 
ubiquitous size distribution found by Zebker et al.l Jl985|) . However, a recent reinvestigation of 
the efficiency of gravitational accretion by means of numerical simulations s hows that this is most 
likely not the case, largely independently of particle collisional properties (IKarjalainen and Salo . 
2004). Because the dispersion velocity of ring particles is so small, even in perturbed regions, this 
implies in turn that, on average, only minute fragments of particles can be collisionally eroded from 
or accreted onto ring particles, leading to very long time-scales for the evolution of the particle size 
distribution, validating the use of Eq. d24l) . at least with respect to this issue. Note that the particle- 
size distribution within the wave may be different from the particle size distribution outside the 
wave, because the coll isional regime is different. In this respect, the spectral "halos" found by 



Nicholson et al.l (I2008|) may indicate a variation in particle sizes in the vicinity of strong density 



waves compared to other ring regions. In any case, the presence of the wave certainly changes the 
particle velocity dispersion and therefore the collisional regime, leading to evolution of the particle 
size distribution on long enough time scales throughout the wave region. This may be accounted 
for by allowing the opacity K = r/S to depend on semi-major axis (any possible dependence on 
the azimuthal direction being quickly erased by the ring velocity shear). 

A related issue concerns the presence of self-gravity wakes, as such structures form and dis- 
solve on time scales comparable to the orbital period (see also IKarjalainen and Said 120041 and 
references therein). It is unclear at this point how the presence of density waves (which locally 
enhance or reduce the density) affect the existence and properties of wakes, and how such wakes 
might affect the optical depth profiles in the wave region. Also, wakes, and more generally regions 
where the particles are closely-packed (as may be the case in sections of the B ring or in density 
wave maxima), certainly affect the way the optical depth is reconstructed from the raw data. 

These last issues are difficult to quantify. We therefore ignore them for the time being, but 
we keep in mind that they may affect in an uncontrolled way the application to real data of the 
reconstruction method developed here. 

In the simulated data presented here, we have assumed for simplicity that the opacity K = 
r /£ = to/£o is constant, so that specifying r D as a function of semi-major axis also specifies the 
radial dependence of the surface density. We also ran simulations with a variable opacity, with no 
impact on the precision of the reconstruction method. 



2.4 Kinematics, and the link between radius and Semi-Major Axis 

We have introduced a number of functions of semi-major axis: e(a), A(a), q(a), 7(a), £ (a), 
To (a). These are the fundamental kinematic functions that our procedure aims at extracting from 
the data; they completely specify the wave kinematics. Note that the recovery of S D from actual 
data requires some dynamical input, e.g., the nonlinear dispersion relation; at the level of preci- 
sion we aim at in the end, uncertainties in the modeling of the ring stress tensor may affect the 
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reconstruction of S , a question we will only briefly touch upon here. 

Furthermore, we have used a and as a set of spatial coordinates; e.g., the surface density of 
the previous section is spatially prescribed in terms of a and <\>. However, the data are given as a 
function of the radius r and azimuth 9. Thus there is a need to apply transformations from a, cf) to 
r, 9 and vice-versa. 

One can check that assuming 9 = in such a change of variable leads to negligible error at 
the level of precision achieved in the data, due to the magnitude of the eccentricity, and the way it 
comes into play in the change of variables (see Appendix lA.il ). Furthermore, in the radio occulta- 
tion profiles, the data acquisition procedure ensures that one can assume that 9 (and consequently 
4>) is constant in a given profiled, with negligible error at the level of precision of the data (see next 
section). Therefore, it is legitimate to assume S(r, 9) = S(r, <p) (or a similar relation for r), but 
one needs to specify how S(a, (f>) is deduced from E(r, <p). 

Eq. (fTTT) allows one to compute the radius r as a function of the semi-major axis a and azimuth 
as long as the functions e(a) and A(a) are known; from our previous remarks, mcf) can be 
considered as a constant parameter depending only on the optical depth profile considered, so that 
this relation can be viewed as an r(a) relation, and can be inverted to find a{r). The details of this 
inversion procedure are discussed in Appendix lA.il 

Note that even though ae is a small quantity compared to the wavelength, it induces a difference 
between r(a) and r(r), which can produce a significant error in the determination of q and r , and 
cannot be neglected. 



3 Simulated Data 

In order for our simulated data to be realistic, we have extracted some of the required parameters 
(such as longitudes and times) from the diffraction-r econstructed occul tation data for the eight 



Cassini Radio Science Team occultation experiments (|Kliore et all |2004|) that will be used later 
on as a real test case of our method. These data were acquired on days 123 (May 3rd), 141 (May 
21st), 177 (June 26th), and 214 (August 2nd), respectively, of 2005. On each day there was both 
an ingress and an egress experiment, for a total of eight profiles. For each profile, the data were 
acquired at the Deep Space Stations (DSS) of the Deep Space Network DSS- 14, DSS-43, DSS-43, 
DSS-63, DSS-14, DSS-14, DSS-63, and DSS-63, respectively. The X-band data (frequency of 8.4 
GHz) are plotted in blue in Fig. [8] with a resolution of 1 km, with one point every 250 meters^, for 
the Mimas 5:3 density wave. This wave, propagating from an inner Lindblad resonance location of 
132301 km from Saturn's c enter, is one of the strongest and most nonlinear waves in the rings. It 



was particularly studied by Longaretti and Borderiesl (11986I) . These data have a free space signal 



to-noise ratic0 of 54 dB-Hz. The instrumental noise is generally very low, although it is larger at 
the peaks. Peaks that reach an optical depth of 5 are not physical because they exceed the RSS 



6 This is not in general true for stellar occultations. 

7 The sampling rate of the diffractio n-limited data is 250 m, and when diffraction is removed the effective spatial 



resolution is 1 km (Maro uf et al . 1986). 



8 Defined as 10 log 10 ( averaged carrier power / average noise power in 1 Hz bandwidth ). 
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capabilities; in other words, we have no data points for the tops of these peaks. In addition to the 
instrumental noise, the data also contain unmodeled physical effects, such as peak shoulders and 
split peaks. 

We consider each occultation of a given feature to be instantaneous with its time t being the 
mean time of the occultation. This assumption introduces negligible error, as can be checked by 
comparing the total data acquisition time for this wave (~ 15 seconds for the Mimas 5:3 density 
wave) with the dynamical time-scale associated with the pattern speed (18.1 hours for the Mimas 
5:3 density wave). We also assume that the ring intercept point of the spacecraft-to-Earth ray 
occurs at a constant mean inertial longitude 6. This assumption is justified for nearly diametric 
occultations, such as the ones considered, which typically cover a longitude variation of 0.01°. 
Therefore, we have one value of rruj) for each occultation. These values are shown in Table [TJ 

We used these values of mcj) and constructed eight simulated occultation profiles in the follow- 
ing way. First we defined the parameters ae(a) and To (a) using the array of radii of the actual 
profiles as semi-major axes and the empirical equations: 



ae(a\ 



T (a) 



0.1 + 2 
0.5 + 0.5 



1 + tanh 
1 + tanh 



a — a c 
a — a. 



a — cln 

Q>\ — CL]\f 



(25) 
(26) 



where a\ (abou|?| 132260 km) is the first (smallest) semi-major axis, and (about 132470) is the 
last (largest) semi-major axis; a c = 132275 km and ( = 50 km. We selected a x = a res as the 
resonance location for the simulated data. These choices are made so that the resulting profiles 
loosely resemble the observed wave. The choice of Tn(a), w hich increases at the beginn ing of the 
wave, is consistent with the calculations of IShu et al.1 (I1985b|) and lBorderies et al.1 (I1986Q . who find 
that the surface density is enhanced in the wave zone outward of the resonance. This is a general 
feature of strong waves. 

Next, we defined the opacity K, and having the background optical depth, this allowed us to 
compute the background surface density S in each point. As mentioned above, the simulated data 
used in this paper have a constant opacity, which is K = 0.0125 cm 2 /g. We also made tests with 
an opacity varying with the semi-major axis. 

The parameters mA and q are next constrained in each point by enforcing the ideal (WKBJ, 
pressureless, dissipation free) nonlinear dispersion relation throughout the wave region. This is 
certainly not satisfied in actual waves, but provides a convenient way to test our ability to recon- 
struct all kinematic parameters, including the surface density, in the wave propagation region. The 
parameters mA and q are computed iteratively in the following way. Before starting the itera- 
tive calculation, we computed in each point gcos7 from Eqs. (fl4l) and (|25l ), and we initialized 
mA(ai) = 0. Next we computed: 



9 This value is in fact the initial radius of our simulated data instead of its semi-major axis. A similar comment 
applies to ajv- 
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fin 

mA(d„) = mA(a n _ 1 ) + / kda, (27) 

J 'a„_i 



where the wavenumber A; comes from the dispersion relation (|Borderies et all 1 1986) 



k = 3(m- 1)Q 2 (a-a res ) _ 
27rGE C(g) a res 

in this equation, f2 is the mean motion at a res , G is the gravitational constant, and: 

s 4 /" +0 ° , sin 2 M r /g 2 sin 2 M\ 



with: 



tftf) = (30) 

Having mA, it is possible to compute g sin 7 from Eq. (fl~5T) and then g. We iterate the calcula- 
tions until the maximum change in the values of q at the same point and between iterations is less 
than 0.01. At the first iteration, C(q) = 1 (i.e., the dispersion relation reduces to the linear one). 

The parameter 7 is computed from Eqs. (fl4l) and (fT5l) : / is then derived from Eq. (f2Th . The sim- 
ulated wave parameters are shown in Fig. [2j Note that 7 increases from to ir/2 for the first thirty 
km of the wave, so that in this region, our simulated data are not consistent with the dynamics of 
a real wave (as the dispersion relation fails there), although self-consistency of the kinematics has 
been enforced in the way the simulated data have been constructed. The simulated data neverthe- 
less resemble an actual wave in the first waveleng th. The failure of the WKBJ ord ering inside the 



first wavelength is consistent with the findings of lLongaretti and Borderiesl (|1986|) for the Mimas 
5:3 density wave. In this portion of the wave and farther, q increases as the wave becomes more 
and more nonlinear until its growth (implied by self-gravitational stresses) is stabilized to a plateau 
and then damped by viscous stresses. Our simulated data also loosely reproduce this behavior. 

Profiles without noise were derived from the above simulated kinematic parameters, the values 
of m0, and Eqs. (fTTI) . (12"2l) . and (l24l) . These profiles are shown by the blue lines in Fig. HI The 
red lines on top of the blue lines represent the reconstructed profiles, which will be discussed in 
section |4~3l 

Next, we applied eight shifts in radii representing our imperfect knowledge of the absolute ra- 
dial scale in actual wave profile recording. These shifts are 0, —0.5, —1.0, 0.5, 1.0, —1.0, —1.5, 1.5 
km. In other words, the first synthetic profile without noise was not shifted, the second one was 
shifted by -0.5 km, the third one by -1.0 km, etc. The ability to recover the values of these radial 
shifts will be one test of our analysis method. 

Next, we generated noisy data from the shifted profiles without noise. For this, we used the 
actual data from 131980 km to 132190 km, i.e. inward the Mimas 5:3 density wave profiles, in a 
region where the optical depth fluctuations are somewhat muted. We could have used free space 
data, but these data contain very little noise. Our assumption was that whatever physical effects 
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give rise to optical depth fluctuations in the region inward of the density wave plausibly acts also 
in the density wave. We used the actual optical depths of each profile in the region considered and 
inferred the amplitudes of the signal. For each profile, we detrended the amplitudes by removing 
a parabola. We computed the standard deviation of the detrended amplitudes. We generated an 
array of pseudo-random values drawn from a normal distribution with a mean of zero and standard 
deviation that is twice the standard deviation of the detrended amplitudes. We averaged over four 
consecutive samples to produce an array of amplitude noise. This technique was used because 
the resolution of the data is 1 km, but they are sampled at 250 m intervals. We added this noise 
to the amplitude of the simulated data without noise, and hence generated noisy simulated data. 
The noise in these data is not i nstrumental noise, it represents a real effect resulting from real 
fluctuations in the optical depth ( Showalter and Nicholson 1990|) . The noisy profiles are shown by 
the blue lines in Fig. [51 Again, the red lines on top of the blue lines represent the reconstructed 
profiles, which will be discussed in section H31 

The simulated profiles were interpolated to a common radial scale extending from 132260 km 
(the assumed resonance radius) to 132470 km, with a step-size of 250 m as in the actual data. 
As mentioned above, our goal was to generate plausible-looking profiles for which the wave pa- 
rameters are known, rather than accurately replicating the Mimas 5:3 profiles. These simulated 
profiles were used to establish the procedure which is described in the next section. We chose 
the example of the Mimas 5:3 density wave to guide us for generating the simulated profiles be- 
cause this wave is clearly visib le and isolated i n the data, and has previously b een studied by 
Longaretti and Borderiesl ( 1986|) . The analysis of Longaretti and Borderiesl (|l986|) used only one 



optical depth profile (a photopolarimeter profile), which did not allow us to accurately constrain 
the location of the peaks, the background optical depth, and the nonlinearity parameter q. In this 
paper, we will apply the procedure to eight actual Mimas 5:3 radio optical depth profiles. 



4 Procedure 

4.1 Outline of the inversion procedure 

Observations (or simulated data) provide us with r(r, 6) = r(r, <p) at specific azimuths 4> ( m the 
frame rotating with the pattern speed). These observations are modeled using r(a, <p) given by 
Eq. (1241) . using the a(r) relation discussed in section [2~4l and Appendix IA. 1 1 so that r(r, 0) = 
r(a(r), <p). They involve three "secondary" functions of semi-major axis: the unperturbed optical 
depth r , the nonlinearity parameter q and the phase function / (Eq. |2D. Furthermore, the quanti- 
ties q and / depend on the "primary" kinematic quantities, i.e. the eccentricity e and the lag angle 
A (one can choose 7 instead), through Eqs. (fl4"l) . (fTBl) . and (|2T|) . Our objective is to recover the five 
functions r (a), q(a), mA(a), ae(a), and 7(a) from the data. The major a priori difficulty lies in 
the fact that the relation between r(r) and r(a) at any given azimuth involves the knowledge of e 
and 7, which are not known beforehand when fitting Eq. ([241) to the data. This makes the recovery 
of a common set of kinematic parameters for all profiles a challenging task. 

One possible method of inversion that would circumvent this problem would be to perform a 
least-squares fit for each radius (and at all radii at the same time), to the five previous parameters 



18 



with the help of the r(a, <p) and a(r) relations. As one has more profiles than unknowns, this proce- 
dure would seem to be well-posed, and would have a definite advantage: no dynamical assumption 
of any kind (such as the WKBJ ordering) would be required, only generically valid kinematic con- 
straints would be used. We have tried this, and, unfortunately, this fails for several reasons. First, 
the number of profiles is limited, and some of them have similar values of mcp, so that they barely 
provide enough independent information. For a simulation, we could have used more profiles and 
a good sampling of values of m0, but we wanted to see what we could do with the Radio Science 
data that have been pre-processed so far. Secondly, the inversion procedure is always extremely 
sensitive to errors in the determination of the phase function /, so that, unless provided with a 
guess unrealistically close to the actual phase function, an iterative global least-squares procedure 
was found to fail systematically. Thirdly, the radial origin of reference of all the profiles is not the 
same, so that one must also fit for the radial shifts one needs to apply to normalize the radial scales. 
In the future we will know the radial absolute scale to within 100 meters, but again we wanted to 
see what we could do with the existing Radio Science profiles. Note that even without radial scale 
uncertainty, the resonance locations vary with time due to the long period variation of the orbits 
of satellites such as e.g. Mimas and Pandora. We found that radial shifts considerably worsen the 
question of the precise determination of the phase function in a global least squares fit. Fourthly, 
the smallness of the wave amplitude and the loss of the WKBJ ordering in the vicinity of the wave 
resonance makes the least-squares procedure less sensitive to variations of mA in this region of the 
wave. Finally, the relations between q and / on the one hand, and e and 7 on the other, involve a 
radial derivative; it is difficult to have these derivatives numerically well-behaved in a global fitting 
procedure. These difficulties make the identification of a common set of kinematic parameters for 
all profiles nearly impossible. Nevertheless, a least-squares global fit may be used as a final step 
to improve the quality of the fit, a point we have not yet checked in any detail due to practical 
implementation difficulties. 

Instead, we have reverted to the approach devised by lLongaretti and Borderies dl986h . with 
a number of significant improvements. The inversion procedure is constructed from the iteration 
of two main stages. In the first one, the secondary quantities r (a), q(a) and f(a) are deduced 
from the data for given primary parameters e(a) and 7(a) (which provide us with a given a{r) 
relationship); this step relies on the WKBJ ordering, which we have shown to be valid at least in 
the nonlinear part of density waves, but we incorporate corrections due to the imperfection of this 
ordering. In a second stage, e(a) and 7(a) are redetermined from the q(a) and f(a) just obtained, 
with the use of Eqs. (fl4l) . (fT5l) , and (I2TT) ; the WKBJ ordering is also used for convenience in this 
stage, but this is not strictly necessary. Determination of the radial shifts of the profile with respect 
to a common radial scale are interweaved in these steps in a manner that will be specified below. 
These two stages are iterated to improve the determination of To (a), q(a) and f(a) on the one hand, 
and e(a) and 7(a) on the other. Iterations proceed until convergence is achieved. We now turn to a 
more explicit and detailed description of these two stages. 

First, an initialization of all quantities is required before any iteration is applied. This initial- 
ization and its rationale are described in the next subsection. 

Next, we present the first stage (recovery of the secondary parameters), which is broken up 
into several steps. The first one (subsection 14.31) describes the determination of the secondary 
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parameters t , q, and ip as a function of the semi-major axis in the vicinity of the peaks for each 
of the various profiles. This is the most complex step of our inversion procedure; it relies on the 
WKBJ ordering. However, since a precise determination of the phase function turns out to be 
critical in our inversion procedure, we apply in this step a correction to the leading order WKBJ 
solution for tp by taking into account the small but non-vanishing variation of the wave amplitude 
compared to the variation of the profile phase; this correction might not always or everywhere 
be useful, but is systematically computed, for simplicity. In the second step (subsection 14.41 ), we 
determine the global phase function / from the profile phase ip as well as a first approximation of 
the radial shift for each profile. We also make use of these profile radial shifts to "unshift" the data 
as well as r , q, and / such that they are on a common radial scale. The last step determines the 
mean secondary parameters r , q and / common to all profiles (subsection !4.5l) . 

The second stage is divided into two steps: recovery of the primary kinematic parameters e and 
7 (and coincidentally, mA) from the mean secondary parameters q and / (subsection !4.6l) . and final 
fitting of the reconstructed profiles r(r, <p) on the data, which allows us to determine an additional 
correction to the radial shift of the profiles (subsection l4.7l) . This correction is needed because first 
we are dealing with an iterative process with several passes through the first and second stages, and 
second, we are making assumptions that are only approximations, e.g., that r and q are constant 
locally around each peak. The wave kinematic parameters t , q, f, ae, and 7 are shifted by the 
same amount as the profiles. 

The WKBJ ordering is primarily needed in the first step of the first stage, and, to the same 
degree of approximation, in the recovery of e and 7. The a(r) relation is needed in the first step of 
the first stage and in the last step of the second stage. Also, as will be shown shortly, this ordering 
implies that using only information at the profile peakfS or in the vicinity of the peaks minimizes 
errors in the knowledge of e and (small) deviations of 7 from ir/2. Because all five kinematic 
parameters we wish to recover vary smoothly enough with semi-major axis, their determination at 
the peaks of the profiles is sufficient to constrain them throughout the wave region. Information 
is therefore gathered in the vicinity of peaks only in the first steps of our procedure. This limited 
information is sufficient to precisely reconstruct the wave kinematics, but naturally, the goodness 
of the fit is checked throughout the wave region by the reconstruction of the profiles. 

In practice, the WKBJ ordering, as already mentioned, is not robust enough in the first wave- 
length region (and in th e evan escent region as well); this is true of actual data, as found by 



Longaretti and Borderiesl (|1986l) . and of our simulated data as well. As will be seen, even in this 
region, our reconstruction method does not entirely fail, but does not provide good enough results 
for our ulterior purpose of torque measurement. Because of this limitation, practical physical diag- 
nostics of wave and disk physics may require us to apply extraneous constraints to this inner wave 
region. This point will be further discussed in our conclusion section. 

Our inversion procedure makes use of three different iterative processes: one for the computa- 
tion of a(r) (see appendix lA.il) . one for the recovery of e and 7 (see appendix IA.2I) . and the entire 
procedure itself. We found that, in general, three iterations of the whole inversion process were 



10 The same argument holds at troughs. However, peak locations are much more precisely defined than troughs, 
although the experimental noise is larger there; only peaks will be used in the procedure outlined in this paper. The 
neglect of the troughs is compensated by the fact that we have many profiles, whose peaks occur at different radii. 
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sufficient to achieve convergence to a high level of precision, with and without noise superimposed 
on our synthetic profiles. Iterations of the whole inversion procedure are referred to as "passes" in 
the remainder of the paper. 

The reader may ask why such an involved inversion procedure is required. We found by trial 
and error that simpler choices either failed entirely, or produced unreliable results. This inversion 
procedure has been implemented in MATLAB. 



4.2 Initialization 

Our initialization is designed such that the kinematic parameters can be treated as the results of a 
fictitious previous pass, for implementation convenience. It assumes that the eccentricity is equal 
to zero everywhere and that the profiles are not shifted. Thus the semi-major axes are equal to 
the radii. As explained below, we focus first on the peak regions, where r = a holds to a high 
level of precision, so initializing e = is a reasonable choice to initiate the procedure. For the 
other parameters, we initially assume 7 = ir/2 (strict validity of the WKBJ ordering), / = (an 
arbitrary choice, but a better one would be required only if e 7^ 0, as it is only used in the a(r) 
relation), r = 1 (from the behavior of r at the beginning of the wave), and q = 0.5 (midway 
between its possible minimum and maximum values). For each profile, and for each peak, the 
peak location is approximated initially by the semi-major axis (or radius) of highest optical depth 
in a window of data containing the peak. 



4.3 Recovery of the secondary parameters: profile phases and peak charac- 
teristics 

We now come to the crucial problem of determining precisely the peak locations as well as the 
corresponding r , q, and ip of each profile using as initial estimates the results of the previous pass 
(or of the initialization in the first pa ss). This is one of the points on which we significantly improve 
on Longaretti and Borderiesl (|l98q) . At each pass through our procedure, we first improve the 
determination of if) (a) (subsection 14.3 .11) , and then the determination of the peaks' characteristics 
(subsection |4.3.2| ). 

In this step, we treat the various profiles independently. 



4.3.1 Profile phase determination: 

Consistently with the WKBJ ordering, we note from Figs.|31H and|5]that the large-scale variations 
of the wave parameters occur over a scale larger than the wavelength. The variations in optical 
depths (see Eq. [24]) are therefore dominated by the variation of ip (see Eq. l22l) with semi-major 
axis. 

Let us first look at what happens if this approximation is exact, i.e., if 7 = n/2 exactly through- 
out the wave region and the variations of q and r with semi-major axis are negligible; these ap- 
proximations cannot hold exactly at the same time (we will come to this below), but constitute a 
useful starting point. In this case, Eq. (1241) implies ip = 2mr at peaks (n is an integer, and successive 
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peaks correspond to successive values of n). It follows that cos (m<p + mA) = cos (tp — 7) = 
and thus r = a at the peaks, so that ip(a) = ip(r) at the peaks. We treat the resonance radius as a 
first "peak" with / = (approximately true at resonance) so that ip = m<p there (modulo 2%). 

It is interesting to note that in the WKBJ approximation, the value of the profile phases at the 
peaks is independent of the peak locations; this is not quite the case if we take into account the inac- 
curacy of the WKBJ approximation. To account for this inaccuracy outside the first wavelength^ 
we compute a correction with the help of the peaks' characteristics determined in the previous pass 
(or of the initialization in the first pass, and the correction vanished then) in the following way. 

Note that because q and r are not exactly constant, the assumption that r(r) = r (a(r)) reaches 
a maximum at points such that r = a is incorrect for several reasons. First, even ignoring the non- 
constancy of t and q, r would indeed be maximum for ip = 2mr, but this would not correspond 
to locations where r = a because 7 = tt/2 would not be exactly true. Secondly, there may be a 
model error in that we have taken r to the lowest order in eccentricity in Eq. (l24l) . Thirdly, r and 
q vary. Even though they vary slowly, these variations imply that the maxima of optical depth do 
not correspond to ip = 2nn. Of all three effects, only the third is possibly non-negligible at the 
level of precision of the data, as can be checked from order of magnitude estimates. The second is 
negligible due to the very small eccentricities (< 1CT 4 , typically). The first is automatically taken 
into account in our procedure in all passes but the first, although it is negligible except close to the 
resonance. The third effect is accounted for by computing a first order correction Sip to the phase 
function, ip = 2nix + 5ip, at the maxima of r(r), due to the non-vanishing derivatives Tq = dr Q /da 
and q' = dq/da. To linear order in 5ip, the constraint dr/dr = at the peaks semi-major axis a pk 
yields: 

x = T o(a P k) [1 - q(a P k)} + r (a pk )q'(a pk ) 
pk q{a P kW{a P k)ro{a pk ) 
This correction is appropriate except close to the resonance where 7 ~ 0. 



4.3.2 Improved determinations of the peak locations, and of q and r at the peaks: 

The procedure outlined in the previous subsection defines ip(a p k) at the peaks of each profile, 
and ip(a) is obtained from ip(a pk ) by interpolation; this procedure is quite accurate as ip(a) varies 
smoothly enough. Next, we make use of this information to improve the determination of the 
positions of the peaks, and to determine the values of q and t at the peaks. To achieve this, we 
rely again on the fact that these two quantities should be approximately constant in the vicinity of 
any peak (or of any point for that matter) and find them by least-squares fitting of the optical depth 
profile in the vicinity of each peak with a function inspired from Eq. (|24j) . 

We first define an approximate optical depth profile function that is valid around any peak a pk 

by 

T f*t(a) = ° pk , ni , (32) 

[1 - q pk cos^(a)J 

"inside the first wavelength, as the WKBJ ordering fails anyway, an altogether different reconstruction procedure 
is required. See the discussion section. 
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where r 0pk and q pk are the values of r and q at the considered peak, and ip(a) is the profile phase 
function just determined. The use of an accurate approximation to the true ip(a) turns out to be 
crucial for the convergence of our inversion method. We then use this approximate representation 
of the exact peak profile while allowing for a possible error 8 pk in the peak location by defining: 

Tfu{r) = T fit (a(r + 5 pk )) , (33) 

which amounts to a radial translation of r(r). This expression involves the knowledge of the func- 
tion a(r), which is obtained from the inversion of r(a) through the process described in Appendix 
lA.ll fin the first pass, e = so that r = a). This inversion uses the determination of e and mA 
obtained in the previous pass. 

This function is then fitted by least-squares to the actual peak profile in a window around the 
considered peak. This fit yields the three unknown parameters of the function: q pk , r 0pk and 5 pk . By 
trial and error, we found that the window was optimal when allowing ip to vary by ±2%/ '3 around 
its peak value; this results from a compromise between the validity of the assumption of constant 
q and t around the peak, and the precision of the method. 

When e ^ 0, T/j t (r) is narrower in the peak regions than Tf it (a) which is a characteristic of the 
actual peak profiles. Due to this property, our procedure results in substantially improved values 
of q pk and r 0pk in the second and subsequent passes with respect to the first. Note that 5 pk is a peak 
radial position correction pertaining to each peak, and is not related to the radial shifts applied in 
section [3] to mimic present inconsistencies in the ring radial scale between actual profile data. 

Finally, we update the peaks' semi-major axes to account for the S pk corrections with the help of 
the a(r) function, so that, for each peak, we now have better approximations of a pk , ip(a pk ), q(a pk ) = 
Qpk, T (a pk ) = r 0pk . The full functions ip(a), q(a) and r (a) are specified by interpolation. 

4.4 Recovery of the phase function and first radial scale correction 

At this point in our procedure, we have determined the profile phase ip(a) for each profile. What 
we need is a single phase function f(a), as precisely determined as possible due to the sensitivity 
of the procedure to errors in this quantity. 

First we compute the phase function for each profile by subtracting rti(f) from ip [see Eq. (1221)1. 
Next we shift the phases by an integer number of 2n. Indeed, the way we have ascribed values 
to ip at each peak does not ensure coherence between profiles since peaks in different profiles 
belonging to the same wavecrest may not correspond to the same value of ip, which may differ 
by an integer number of 2n. Once this is done, the various / functions still present some scatter, 
which is due both to errors in our reconstruction procedure and to differences in the absolute radial 
scales of the various profiles. We determine a first approximation of these absolute radial scale 
differences by translating each profile by an integer multiple of 2tc to minimize the difference 
between the phases of the current translated profile and the first profile. In this way we compute 
a single fixed radial shift for each profile. We call this action "collapsing the phases" because if 
there were no approximation or assumption in our whole procedure, all the phase functions would 
coincide after this collapse. This determination is of course relative, since the various radial shifts 
are computed by arbitrarily taking the first profile as a reference profile. To compute the radial 
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shifts, we limit ourselves to the range from a = 132320 km to a = 132440 km, where the phase 
is well-behaved and little affected by boundary effects. The phases agree to a few percent in the 
range from a = 132320 km to a = 132440 km, which we consider a very good result. 

The radial shifts (or more precisely "unshifts") in the semi-major axis (or radius) just deter- 
mined must be accompanied by similar unshifts of r and q, f, and of the data themselves. This 
requires some interpolation over an interval in semi-major axis between 132270 km and 132440 
km. 

4.5 Recovery of the secondary parameters: mean solution construction 

We are now in position to use the functions r (a), q(a), and f(a) for each of the eight profiles 
to produce a single function for each quantity that should be valid for all eight profiles. In other 
words, we want to determine the mean solution for the functions r , q, and /. We treat each of 
these three functions separately. 

It turns out that a polynomial fit to a simple average of each function on each point does not 
satisfactorily represent the true function (i.e., the function used to generate the simulated profiles). 
Instead, we found that the following least square procedure yields good results. The variables of 
this fit are the mean values of the function over a limited number of points or nodes (we chose 
10 for this wave) approximately regularly spaced between 132270 km and 132440 km, a range 
of values determined empirically to avoid edge effects. These variables are first estimated as the 
average values of the eight functions at each node, providing initial conditions to the least square 
fit. The function we wish to fit is computed everywhere by interpolation over these 10 values, 
and is constrained to best represent, in the least square sense, the set of eight profile-dependent 
functions between 132270 km and 132440 km. In order to avoid negative, unrealistic values of 
q, this quantity is constrained to a small, positive value at the resonance radius. This somewhat 
artificially constrains the behavior of q in the first wavelength. The negative values of q derived 
from the various profiles arise in part from the failure of the WKBJ approximation. 

Figure [3] shows the mean solution in red for r , q, and (/ — f injected) /(2tt), on the six upper 
panels. The left panels correspond to the case of data without noise, while the right panels cor- 
respond to noisy data. The green lines refer to the injected simulation parameters (i.e. the wave 
parameters used to generate the simulated data and shown in Fig. [2]>. The blue lines show the 
values determined independently for each profile. All the plots in the six upper panels are cut at 
132270 and 132440 km because the failing of the WKJB approximation at the beginning of the 
wave and boundary effects near the end of the data set make the solution unreliable outside these 
cut-off values. As a matter of fact, the solution is incorrect within the first 50 km of the wave, 
where the wave function / is not determined accurately. As expected, there is more scatter in the 
case with noise than in the case without noise. Nevertheless, the agreement between the mean 
solution and the injected values between 1323 10 and 132440 km is remarkable. This illustrates the 
efficacy of this approach to find the mean solution for r , q, and /. The four lower panels will be 
discussed in the next section. 

The improvement that this method brings to the determination of r and q, as compared to the 
individual determinations, can be quantified by looking at Fig. [6l The left panels are for the case 
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without noise and the right panels for the noisy data. The top panel of this figure pertains to r 
and the bottom panel pertains to q. The middle panels will be discussed in section |4~8l In each 
panel, the magenta lines display the function used to generate the simulated profiles (or injected 
function). The green and cyan lines show the minimum and maximum values F min and F max , 
respectively, of the computed F = r or F = q at any semi-major axis for the eight profiles. The 
red lines represent the mean solutions computed in this section. The black areas represent one 
standard deviation error bars computed as (the "noise" is statistically independent in the various 
profiles): 

{Fmax F m i n ) 

° F = 9 m — » (34) 

MpROF 

where Npro f = 8 is the number of profiles. We see that for both r and q, the magenta lines, which 
represent the true solutions, fall within the error bar areas, and most of the time coincide with the 
mean solutions, shown in red. Therefore Eq. (l34l) provides a reliable handle of the error on tq and 
q. Also this figure demonstrates the need to have as many profiles as possible (preferably well 
distributed in m0) to compute accurately r and q. A precise determination of these parameters is 
especially important to determine the ring surface density and to constrain the ring stress tensor, 
which will be done in a future paper. 



4.6 Recovery of the primary parameters: determination of ae and 7 

To solve for ae and 7, we solve iteratively Eqs. (fl4l) and (fl~5l) . a procedure described in appendix 
IA.2I The equations are solved between semi-major axes of 132290 km and 132440 km, as this iter- 
ation procedure involved in the resolution relies on the WKBJ approximation, which is not robust 
enough in the inner wave region. Where needed in other subsequent passes, ae and 7 are linearly 
extrapolated outside the range specified above. This extrapolation is certainly one major source 
of error in this reconstruction procedure, which explains the relatively large (albeit tamed) differ- 
ence between the injected and reconstructed phase functions within the first wavelength region. 
Improvement is clearly required on this point. This is discussed further in our last section. 

Figure [3] shows the mean solution in red for ae and 7 — n/2, on the four lower panels. As 
mentioned above, the left panels correspond to the case of data without noise, while the right 
panels correspond to noisy data. The green lines refer to the injected simulation parameters. We 
find a remarkable agreement between the mean solution and the injected values over most of the 
wave. There is a difference in the first 50 or 60 km of the wave, where the phase function was not 
well determined due to the failure of the WKBJ approximation. 



4.7 Profiles reconstruction, and second correction to the radial scale 

We are now in position to construct optical depth profiles r(r, <p) = r(a(r), <p) over the entire radial 
range of the eight profiles, to check our obtained kinematic solution. We use Eqs. (fTTI) and ([24]) . 
In doing so, we determine by least squares a second approximation to the radial shifts required by 
the inconsistent absolute radial scale of the various profiles and the various approximations made 
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in the procedure. As before, this second shift translates into a correlative shift of all our kinematic 
functions. The final radial shifts are shown in Tables CD The differences between the applied and 
the computed data shifts are in general significantly smaller in the absence of noise but even with 
noisy data, they amount to at most several tens of meters, which is much smaller than 250 m, the 
separation between the data points. 

Figures |4] and [5] show the unshifted reconstructed profiles (in red) together with the unshifted 
simulated data (in blue) for the profiles without (left panels) and with noise (right panels). We can 
see that the agreement is excellent, and suggests that a better determination of the radial scale shifts 
is hardly possible. Note that the fit is still tamed inside the first wavelength region, although the 
reconstruction of ae, 7, and mA is only correct within a factor of two in this region. Nevertheless, 
an improved determination of the kinematic parameters will be required in this region to compute 
satellite torques with some accuracy. The very last peaks are not accurately represented because the 
fit was not performed in this region; the reconstructed parameters have simply been extrapolated 
there. 



4.8 Recovery of the surface density 

Recovering the surface density requires some further dynamical input, which is usually provided 
by the dispersion relation. As the effect of the ring stress tensor is small, the pressure correction is 
usually neglected in front of the self-gravity term in the dispersion relation, which then reduces to 
the following form in the region of validity of the WKBJ ordering: 

K _ 2irGC(q)kT a res 

3(m — l)Q 2 (a — a res ) ' 

where K = r /E is the opacity. This expression wa s applied in the construc tion of our simulated 
data. The situation with real waves is more complex. Sremcevic et al.1 l2008i) have argued that the 



pressure term leads to measurable changes in the dispersion on the order of 5%. For our application 
to the Mimas 5:3 density wave, we have taken this term into account [see Eq. (|39|)1. 

The opacity is the quantity we chose to recover, as the optical depth r , the nonlinearity pa- 
rameter, and the wavenumber k are known from the inversion procedure; k is computed as the 
derivative with respect to the semi-major axis of the phase function /. 

We also make use of the error estimates on r Q and q to estimate the error on the recovery of the 
opacity K and the surface density S . In doing so, we can neglect the error in k because the error 
in the slope of / is smaller than the error in /, leading to relative errors on k of the order of one or 
two percent at most, much smaller than the error in the two other quantities (of the order of 10 to 
20%). 

We have found that the deviations in q and r a of any given profile with respect to the mean are 
not independent, but rather tend to be anti-correlated. Therefore, they were not added. Instead, we 
computed K for five different cases, involving several combinations of To and q. Remember (see 
section H3T) that we denote by To >m i n (or q m in) and r 0:max (or q ma x) the minimum and maximum 
values of r (or q) in each data point computed independently for each profile. We denote by r 
and q the mean values of r and q computed in section [431 
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The values of the opacity are shown in the medium panels of Fig. |6j The left side corre- 
sponds to the data without noise, and the right side to the noisy data. The nominal value of the 
opacity, K(r ,q) is displayed in red while the injected value of K = 0.0125 is plotted in ma- 
genta. In the second panels from the top, the green lines show K(ro )m i n ,q) and the cyan lines 
show K(T 0tmax ,q). In the third panels from the top, the green lines show K(t , q m %n) and the cyan 
lines show K(r , q ma x)- The error bars, drawn in black, are computed by an equation similar to 
Eq. (|34l) . One can see that the true solution (magenta), stays in the error bar area of the nominal so- 
lution between 60 and 170 km of the beginning of the wave. This demonstrates that the procedure 
leads to a much better solution than that could be inferred by considering the green and cyan lines, 
and that the black area (one standard deviation) gives an accurate estimate of the reconstruction 
error of the opacity. 

The surface density can now be recovered easily as S = tq/K. It is plotted with red lines on 
Fig. U\ The top panel is for the case without noise and the bottom panel is for the case of the noisy 
data. The green lines represents the injected value. 

We compared our background surface density with what would be predicted by the linear theory 
for each profile independently. For each profile, we determined k as the derivative of the profile 
phase ip [see Eq. ((22)) 1. and we computed a constant background surface density by a least-squares 
fit of the equation: 

a - a r es _ 2nGk 
a res S 3(m — 1)Q 2 ' 

which is the linear dispersion relation. This is similar to what is routinely done in the literature. 

The resulting values of S for each profile are displayed as blue lines on the two panels. These 
values are above the average of the values obtained by using the nonlinear theory with our proce- 
dure. Fig. |7] clearly demonstrates the crudeness of the linear theory when the surface density is 
not constant, which is expected in real waves (see section [3]). Our ultimate goal requires that we 
recover the variation of the surface density in the wave region; this makes the use of a procedure 
of reconstruction of the wave kinematics such as the one presented here a necessary first step. 



5 Application to the Mimas 5:3 density wave 

The application of our procedure to the actual data collected for the Mimas 5:3 density wave 
raises a number of other difficulties. The first one is related to the variations in the orbital motion 
of Mimas. Table |2] lists the four resonance locations that can be associated with the position of 
Mimas during the four occultations which provide the eight profiles we have analyzed. For each 
date, we have fitted over a period of two days the epicyclic elements to the orbital data given by the 
SPICE kernel sat267, and then we have solved for the position of the resonance. We see that the 
resonance location varies by 2.6 km, i.e. a sizeable fraction of the (varying) wavelength, leading to 
potentially non-negligible modulations of the wave between various profiles. Two limiting cases 
can be considered: in the first, the wave reacts in a quasi-static way to the variation of the orbit of 
Mimas; in the other, the satellite orbital variations are fast compared to the time-scale of adjustment 
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of the wave. Orders of magnitude estimates of the density wave group velocity show that it is 
larger than the sp eed of varia tion of the resonance location, but only by a factor of a few: the group 



velocity ttG^/k (|Shulll984f) is approximately 20 or 30 km/year for the Mimas 5:3 density wave, 
while the semi-major axis of Mimas varies by 2 or 3 km with a period of 225 days (in addition to 
the ~ 70 year period due to the resonance with Tethys). Nevertheless, we deal with this difficulty 
by setting the resonance radius at 132301 km and assuming one can absorb any modulation of 
the wave in the radial shifts of the profiles, as our aim is more to illustrate the behavior of our 
reconstruction method than to derive precise quantitative results. We found empirically that the 
radially rescaled profiles obtained in this way are reasonably consistent with one another, indicating 
that any spatial and temporal modula tion of the wave behavio r remains mild enough, a conclusion 



consistent with the recent findings of IStewart and Sremcevic (2008). 



The second difficulty, already mentioned in the introduction, is t hat self-gravityl wakes cause 



the ap parent background optical depth to vary with viewing geometry (Colwell et al.1120 06: Hed man et al 
2007). We normalized the optical depths of the eight profiles as follows. First we determined the 
mean optical depth Tj for each profile in an apparently unperturbed region of 20 km inward of the 
Mimas 5:3 density wave (see Table [2]); it appears that modulations in the mean optical depth are 
rather moderate. Thus, we computed the average r ave of the Tj and we multiplied the optical depth 
of each profile by r ave /Tj. Nevertheless, the justification of this procedure is unclear as we do not 
know if wakes inside the wave, if present, affect the optical depth in the same way as outside it. In 
any case, this is not a large effect, probably less important than the uncertainty related to the large 
intrinsic noise observed in the radio data. 

Next, we applied our procedure to the optical depth data. Fig. [8] shows the data (in blue) and the 
reconstructed profiles (in red) after three passes. The agreement between the model and the data is 
very good, considering the large fluctuations of the optical depth (notably at the peaks) present in 
the data. Fig.[9]is the equivalent of Fig. [3] with a few differences. In the absence of injected value, 
the phase function is compared to the its mean solution. The dotted lines correspond to the region 
where the functions for at least one of the profiles were extrapolated. In this extrapolation region, 
the functions t and q increase. These functions are nearly constant between 50 and 120 km from 
the beginning of the wave (at 132301 km, the selected resonance radius). In this domain (132351 
to 132421 km), the phase function / is very well determined. As expected, ae increases and then 
decreases to zero, and 7 is very close to n/2. Table [2] lists the shifts found for the eight profiles 
(they are all positive, but this is only a coincidence). 

Fig- [TO] is the equivalent of Fig. [6] but does not show the minimum and maximum values that 
were used to compute the error bars, since they do not provide additional information. The plots 
were traced in the domain where we can most trust the solution, i.e., between 132350 km and 
132420 km. As we have not tried to determine a position-dependent opacity from the data, the 
opacity is not shown in this plot. 

Instead, the calculation of the opacity along with the velocity dispersion is illustrated by Fig.lTTl 
This figure was constructed as follows. We wrote the nonlinear dispersion relation in the form: 

9th = 9obs, (37) 

with: 
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2nGaC(q)kr 

9obs ~ 3(m - 1)^ ' (38) 



9th = ctK 



[a — a r 



2k 2 H(q 2 )c 2 
3(m- l)fi 2 _ 



(39) 



where a = 1 x 10 12 is an ad hoc numerical factor designed to bring g obs and to values of 
order unity and a res is the resonance semi-major axis. The second term in Eq. (1391 is the pressure 
term, in the hydrodynamic approximation where the (vertically integrated) pressure is given by 
a simple isothermal equation of state p = Ec 2 ,; it introduces the isothermal sound speed (of the 
order of the particles' velocity dispersion) c as a second parameter besides the opacity in the 
dispersion relation. Indeed, we found that the data were good enough to allow us to account for the 
pressure correction to the dispersion relation, which leads to a better agreement between g th and 
g Q bs- At this point, it is unclear if the remaining deviations from the purely self-gravitational term 
are due to imperfection in the radio data that would intrinsically bias the determination of q and 
To (a point further discussed below), imperfections in the model of the pressure term, neglect of 
the s atellite term, or position depend ence of the opacity and velocity dispersion; note in any case 



that Longaretti and Borderiesl (|l986l) pointed out that a constant opacity provided a better fit to the 
dispersion relation than a constant surface density. 

We chose to determine only the opacity and velocity dispersion, and not the resonance radius, 
for reasons also discussed below. We find the velocity dispersion c and the opacity K by a least- 
squares fit of Eq. (|37l ) over the interval between 50 and 120 km from the beginning of the wave, in 
the following way. 

In Fig.QTl the solid lines represent the values of: 

9l = aK ^ a - a ^\ (40) 

& res 

and the dashed lines represent the values of: 

_ 2airGT kC(q) 2ak 2 H(q 2 )c 2 K 
92 ~ 3{m-l)Q 2 3{m-l)Q 2 ' ( } 

The black lines were obtained by using the mean values of r and q. The red and blue lines on 
the top panel were obtained by using the mean value of r and the minimum and maximum values 
of q provided by the error bars on this variable. On the other hand the red and blue lines on the 
bottom panel were obtained by using the mean value of q and the minimum and maximum values 
of t provided by the error bars on this variable (one standard deviation). These error bars indicate 
that deviations from strict linearity may be real; this point will be further investigated elsewhere. 
In the nominal case, the solution gives c = 0.55 cm/s and K = 0.018 cm 2 /g. 

Because of the non-constant background optical depth, the constant opacity results in a non- 
constant surface density, shown in red in Fig. [121 The blue lines were obtained for the use of 
the dispersion relation in the linear limit (q <C 1). Note that we have C(g)S ,jvL £o,l where the 
subscripts NL and L refer to nonlinear and linear, respectively. Since C(q) > 1, we expect the 
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nonlinear background surface density to be smaller than the linear background surface density, as 
observed. 

It is important to note that the cutoff of the radio data at an optical depth of about 5 for the 
profiles considered, due to the instrumental noise in the radio data, flattens the peaks to this value. 
This affects the determination of r and q. Also, the unexpectedly large fluctations in optical depth 
observed at the wave peaks with respect to troughs lead to a similar systematic effect. Similarly, the 
noise level precludes a reliable determination of the weak peaks in the wave damping region, which 
is eventually needed to constrain the ring stress tensor. As mentioned earlier, we have verified that 
the errors in r and q are anti-correlated. Attempts to fit Eq. (1371) with other combinations of r and 
q led to the result that the best estimate of a res was obtained for the minimum value of r and the 
maximum value of q. This indicated that r is overestimated and q is underestimated, and led us to 
enforce the position of the resonance in the fit. We conclude that the variations of the background 
optical depth and surface density and the variations of the nonlinearity parameter derived from the 
radio data are imperfectly determined from the radio data. A better determination, which could 
come from using optical occultation data, is needed to study the damping of the wave and extract 
more (and more reliable) information from the dispersion relation. 



6 Conclusion 

We have established a robust procedure to analyze nonlinear density waves and we have applied 
it to the Mimas 5:3 density wave. In principle, the number of independent profiles provided by 
the Cassini mission would seem to make it possible to devise an inversion procedure relying only 
on the kinematic characteristics of the wave, but we found that in practice, and relying only on 
the radio data, this fails for three major reasons: the differences in the absolute radial scales of 
different profiles, the complex way in which the a{r) relation comes into play, and the sensitivity 
to errors in the dependence of the wave phase on semi-major axis. 

Ins tead, we have extended and considerably improved the method introduced in lLongaretti and Borderies 



(|1986|) . which relies not only on the wave kinematics, but on the WKBJ ordering as well. The in- 
trinsic nonlinearity of the problem makes the procedure complex. Nevertheless, it is possible to 
reconstruct the wave parameters as well as the data, in the presence of uncertainties in the absolute 
scale of about ±2 km and of noise in the data; the data resolution adopted was 250 m but the results 
should not be sensitive to this parameter. 

As such, this inversion procedure performs very well in the far wave region. This opens new 
possibilities of physical diagnostics of wave regions. In particular, with the help of the nonlinear 
dispersion relation and the wave damping equation, one may possibly constrain the surface density 
and ring's stress tensor in a joint way, which would provide us with an indirect window on the ring' 
collisional physics. We plan to make use of this inversion procedure to undertake such systematic 
studies in the future. 

The major inaccuracy of this reconstruction procedure is its limited efficacy within the first 
wavelength of propagation of the wave, where the WKBJ ordering is expected to fail on theoretical 
grounds, an expectation confirmed by the behavior observed in actual wave profiles. This limitation 
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is most drastic in the reconstruction of the eccentricity e(a) and lag angle A(a). This is unfortunate 
for the purpose of the determination of the exchange of angular momentum with the satellite, as 
the generic expressio n of the torque results from a r adial integral involving these two quantities 
[see, e.g., Eq. (35) of LLongaretti and Borderieslll986ll . The results of the present inversion method 
(which yields a fairly educated guess of the form of the kinematic parameters even inside the 
first wavelength, and has allowed us to remove the problem of inconsistencies of the absolute 
radial scale between the various profiles), may ensure the convergence of a refined, direct inversion 
method of the type we initially tried, provided enough high quality profiles are available. This 
approach may require to a different method of inversion of Eqs. (fl4l) and (TTBI) . if needed; a possible 
way to do this is to re-express these equations as a set of nonlinear algebraic equations, and solve 
them by an algebraic iteration method. Additionally, further dynamical constraints may be needed 
to ensure that such a direct inversion method is well-defined. Three levels of such constraints are 
available in the literature. In their most useful fo rm, they are all related to the dynamical equation 
governing the behavior of Z = e exp (zrnA) (see lShu et al.ll 1985c ): they are the dynamical equation 
itself, its solution in the linear limit, and the asymptotic behavior of this solution in the evanescent 
region. Defining the most appropriate strategy on this issue requires substantial additional work, 
which will be reported elsewhere. 

The application to the actual data of the Mimas 5:3 density wave gives good results, although 
the determination of r and q is impaired by the noise in the radio data; better results are expected 
from the use of, e.g., the UVIS data. 



A Appendix 

A.l Radius versus semi-major axis inversion procedure: 

To compute the semi-major axis a a function of the radius we make use of Eq. (fTTI) in an iterative 
procedure. The (n + l)th approximation of a, denoted a^- n+1 \ is related to the nth one through 

fl (n+l) = r + a (n)g( a (n)) cog [ m + m A(a (n) )]. (42) 

Convergence is ensured because e C 1; the procedure is initialized with = r. 

The reader may ask why Eq. (fTTT) cannot be analytically inverted, at least approximately. This 
would be a definite advantage in the numerical implementation of our procedure, as the previous 
inversion is quite heavily required, and slows down the whole process. In principle, one could do 
this through an expansion in eccentricity. The zeroth order solution is: 

a = r + re(r) cos(m0 + mA(r)). (43) 
We look for the first order solution by writing a = r + 8\ with 8\ = 0(e). We find: 

re(r) cos(md) + mA(r)) 

^ = — i — hr k\ ■ (44) 

1 — q(r) cosip(r) 
The second order solution is of the form a = r + 5i + 5 2 with: 
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1 [re(r) cos(m</> + mA(r))] 

2 [1 — q{r) cosip(r)] 3 



2 



dq{r) dip(r) . 

— : — cos w(r) — q — ; — sin ip[r ) 
dr dr 



(45) 



Figure [13] shows the error between the approximate a{r) obtained when using the zeroth, first, 
and second order solutions and the actual a(r). Naturally, the actual functions e(a) and A(a) 
(evaluated at a = r) are used in the process. Although the eccentricity is very small, the error is 
substantial, of order 1 km, and the convergence is very slow due to the various powers of the 1/ J 
factors involved in this expansion, and more crucially due to the fact that the order of the A;th order 
term in the expansion is not of magnitude e k but ae(ae/^) k ^ 1 , where £ is the scale of variation of 
the wave amplitude introduced in section [2T2l the factor ae/£ is substantially larger than e. This is 
why we adopted the iterative method described above. 

However, using the same type of eccentricity expansion, one can check the statement made 
earlier about the impact of the 6 = <p assumption; the difference lies in the fact that in this case the 
expansion at order k is really of magnitude e k . 

A.2 Recovery of the primary kinematic paramaters 

The approximate WKBJ ordering which applies over nearly all the wave propagation region sug- 
gests that Eqs. (fl4l) and (IT3T ) can be used in an iterative way to determine e and 7. In practice, 
defining the nth order approximation of e and 7 as and 7^, we use: 

ae (n) = g 5 " 1 ^ (46) 

[df/da-d^/da] 7 ( } 

and 

(„,-,) ldae^ 

cos7 (n+1) = — . (47) 

q da 

Note that these equations are used alternatively and not simultaneously: gives from 
Eq. (|46l) . and in turn gives ^ n+l ^ from Eq. (|47l) . The iteration is started with 7^ = n/2. 
Obviously, this fails inside the first wavelength, as convergence is ensured only when 7 is not 
substantially different from w/2. 

These expressions involve derivatives of / and 7 with respect to a. We compute these deriva- 
tives by using a degree 5 polynomial interpolation of / and a simple procedure for the derivatives 
of ae and 7. Furthermore, we must smooth ae and 7 after each iteration in order to avoid propa- 
gation of numerical noise through successive numerical derivatives. We stop the iteration process 
when the results become stable. 
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Profile # 


1 


2 


3 


4 


5 


6 


7 


8 


171(f) 


-29.477 


-160.081 


176.166 


44.460 


122.805 


16.740 


60.007 


-9.538 


Applied shifts (km) 





-0.5 


-1.0 


0.5 


1.0 


-1.0 


-1.5 


1.5 


Computed shifts (km) 





-0.500 


-1.001 


0.502 


1.002 


-0.980 


-1.441 


1.501 


Differences (m) 








1 


2 


2 


20 


59 


1 


Applied shifts (km) 





-0.5 


-1.0 


0.5 


1.0 


-1.0 


-1.5 


1.5 


Computed shifts (km) 





-0.519 


-1.038 


0.460 


0.978 


-1.009 


-1.524 


1.424 


Differences (m) 





19 


38 


40 


22 


9 


24 


76 



Table 1: First line: Values of mcj) (m = 4) in degrees for each of the eight profiles. Next three 
lines: Computed shifts of the simulated profiles without noise compared to the shifts applied in 
generating the simulated data without noise; absolute differences between the applied shifts and 
the computed shifts. Last three lines: Computed shifts of the noisy simulated profiles compared to 
the shifts applied in generating the simulated data; absolute differences between the applied shifts 
and the computed shifts. 
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Profile 


Rev. 


Date 


CL res 


r 


Radial shift 


1 


7 


May 3, 2005 


132301.717 


0.84 





2 


7 


May 3, 2005 


132301.717 


0.72 


1.550 


3 


8 


May 21,2005 


132302.037 


0.78 


0.564 


4 


8 


May 21,2005 


132302.037 


0.77 


1.449 


5 


10 


June 26, 2005 


132300.601 


0.72 


0.361 


6 


10 


June 26, 2005 


132300.601 


0.77 


0.146 


7 


12 


August 2, 2005 


132299.480 


0.72 


0.830 


8 


12 


August 2, 2005 


132299.480 


0.75 


0.177 



Table 2: Location a res of the Mimas 5:3 resonance for the four dates of each profile. Each date 
has an ingress and an egress so we have eight profiles. Average optical depth f for each profile in 
an unperturbed region of 20 km inward of the Mimas 5:3 resonance. Radial shifts obtained by the 
procedure. 
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r = a[l — e cos(m0 + mA)] 

(9^0 = 

« = -A 



Figure 1: Streamlines of a 4-lobe (m = 4)density wave, such as Mimas 5:3. The lag angle A provides the 
relative shift in angle of one streamline to the next. The loci of constant profile phase ip corresponding to 
the density maxima (tp = 2irn with n =1 to 5, solid line) and minima (ift = it modulo 2ir, dashed line) are 
also shown to illustrate the geometrical meaning of tp. A streamline has been singled out, to show the phase 
angle A. 
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Figure 2: Kinematic wave parameters for the simulated profiles. These parameters were computed as 
described in section [3] The assumed resonance radius is 132260 km. 
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Figure 3 : Solution for the five wave kinematic parameters after three passes for the simulation case without 
noise (left panels) and with noise (right panels). The blue lines show the solutions for the profiles taken 
independently from each other. The red lines display the mean solution. The green lines refer to the values 
injected in the simulation (see Fig. 
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Radius -132260 km 
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Figure 4: Simulated profiles in blue for the case without noise, and reconstructed profiles in red. For each 
profile, the optical depth scale varies between and 4. Intermediate tick marks represent both the r = 
level of the next profile, and the r = 4 level of the previous one. 
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Radius -132260 km 

Figure 5: Simulated profiles in blue for the case with noise, and reconstructed profiles in red. For each 
profile, the optical depth scale varies between and 4. Intermediate tick marks represent both the r = 
level of the next profile, and the r = 4 level of the previous one. 



42 



60 80 100 120 140 160 180 



60 80 100 120 140 160 180 



™ E 0.014 

& 0.012 

o 
cd 

Q. 

o 



0.01 



N E 0.014 
0.012 

o 
cd 

Q. 

o 



60 80 100 120 140 160 180 



0.01 




100 120 140 160 180 



^ E 0,014 

~ 0.012 

o 
cd 

Q. 

o 



0.01 



D5 

N E 0.014 

~ 0.012 

o 

cd 

CL 

o 



60 80 100 120 140 160 180 



0.01 




60 80 100 120 140 160 180 




60 80 100 120 140 160 
Semi-major axis - 132260 km 



180 




180 



Semi-major axis - 132260 km 



Figure 6: Background optical depth, opacity, and nonlinearity parameters for the simulated cases without 
noise (left panels) and the cases with noise (right panels). In each panel, the red line represents the solution, 
surrounded by its error bars in black. The injected values used to generate the simulated data are displayed 
by magenta lines. On the top and bottom panels, the green and cyan lines show the minimum and maximum 
values of the function obtained by analyzing the profiles independently of each other up to the point where a 
mean solution is determined. There are two panels for the opacity corresponding to two ways of computing 
the error bars (see section l4~8l) . This section also explains how the green and cyan lines were obtained. 
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Figure 7: For the simulation, surface density computed by a fit to the linear theory (blue lines), computed 
by the nonlinear theory (red), and injected in the simulation (green). The top panel is for the case without 
noise, and the bottom panel is for the case with noise. 
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Figure 8: Radio optical depth profiles for eight occultations (in blue) of the Mimas 5:3 density wave and 
solution after three passes (in red). For each profile, the vertical scale goes from to 4. Intermediate tick 
marks represent both the r = level of the next profile, and the r = 4 level of the previous one. The 
resonance radius is 132301 km. 
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Figure 9: Retrieved kinematic parameters of the Mimas 5:3 density wave. The kinematic parameters 
obtained by considering the profiles independently from each other are shown in blue. The mean solution 
is displayed in red. The dotted lines refer to the region in which the parameters for at least one profile are 
extrapolated. 
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Figure 10: Mean solution (in red) for the background optical depth and the nonlinearity parameter sur- 
rounded by the error bars in black, for the Mimas 5:3 density wave. 
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Figure 1 1 : Functions g\ and g2 given by Eqs. (l40l) and (l4Tb for various combination of tq and q. Note that 
T~o : min, To,max, Qmin, and q m ax correspond to the values at the boundary of the error bar areas, and not to the 
minimum or maximum values of these variables obtained by treating the profiles independently. 
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Figure 12: Surface density from the nonlinear theory (in red) and from the linear limit for each profile 
considered independently (in blue) for the Mimas density wave. 
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Figure 13: Error between the approximate a(r) obtained when using the zeroth, first, and second order 
solutions and the actual a(r). 
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